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Strong field ionization is difficult to treat theoretically due to the simultaneous need to treat 
bound state dynamics accurately and continuum dynamics efiiciently. We address this problem by 
decomposing the time dependent Schrodinger equation (TDSE) into an inhomogeneous form, in 
which a precomputed bound state acts as a source term for a time dependent tunneling component. 
The resulting theory is equivalent to the full TDSE when exact propagation is used, and reduces 
or eliminates a major source of wavefunction error when propagation is approximated. The gauge 
invariance of the resulting theory is used to clarify an apparent gauge dependence which has long 
been observed in the context of strong field S-matrix theory. 



The ionization of an atom or molecule by an intense 
laser field is one of the basic phenomena of strong field 
physics, and also one of the most difficult to describe 
theoretically. The difficulty arises from the very different 
theoretical techniques needed to accurately describe the 
bound states of the molecule on one hand, and to effi- 
ciently describe the dynamics of the ionized electron on 
the other. In the zero field limit, molecular states can 
be well described using the methods of quantum chem- 
istry. However, such methods typically scale poorly with 
system size, making them ill suited to problems in which 
one electron may travel dozens or hundreds of atomic 
units from the ionic center. Meanwhile, the dynamics 
of the ionized electron can be well described using sin- 
gle active electron calculations, at the cost of giving a 
poor description of the bound states. In the tunneling 
regime, the continuum component is typically exponen- 
tially suppressed with respect to the bound component, 
so that even slight errors in the treatment of the initial 
state have the potential to swamp the desired tunneling 
wavepacket. 

To date, most analytical methods of strong field ion- 
ization are built around a single active electron picture, 
including S-matrix theory^!.] and its variants 0,0], includ- 
ing PPT theory [4] and strong field eikonal-Volkov [l, Q . 
Methods built around quantum chemistry approaches 
include multiconiiguration time dependent Hartree and 
Hartree Fock [THSl and time dependent R-matrix ap- 
proaches 13- U- Coupled channel approaches 14- lij 
allow interactions with the departing electron to drive 
transitions between different states of the ion, so that 
the electron "hole" may occupy any of several orbitals. 

This paper addresses the problem of strong field ion- 
ization by separating the evolution of a time dependent 
"tunneling" component of the wavefunction from that of 
an initial bound state. Because the time propagator does 
not act upon the initial state, it may be approximated 
without fear of erroneously projecting the initial state 
into the continuum. Likewise, the initial state may be 
found using any desired level of theoretical approxima- 
tion without increasing the computational cost of prop- 



agation. The resulting general method reproduces either 
the full time dependent Schrodinger equation (TDSE) 
or strong field S-matrix theory in the appropriate lim- 
its, and is manifestly gauge invariant. By examining 
its transformation under a change of gauges, we clarify 
a longstanding apparent gauge dependence in S-matrix 
theory. 



THE INHOMOGENEOUS SCHRODINGER 
EQUATION 

When subjected to an intense laser field, an electronic 
wavefunction begins to deviate from its field free state. 
If the initial wavefunction is an eigenstate of the field 
free Hamiltonian ijji,{x,t) ~ ipi,{x)e~'^^''* , where Hoipi, = 
Eb^pb, the full wavefunction may be written as the sum of 
this term plus a time dependent "tunneling" component 
(despite the name, the logic works equally well in the 
multiphoton regime) 



(1) 



Mathematically, the tunneling component serves to can- 
cel the nonzero action residual 



(2) 



which results when the field-free eigenstate is acted upon 
by the field-on Hamiltonian. Writing the full TDSE as 



d 



(3) 

and subtracting the evolution of the bound state due to 
the field free Hamiltonian 



-iEot _ 







(4) 



yields an inhomogeneous equation in which the laser field 
H{t) — Hq acts upon the bound state to serve as a source 
term for the tunneling component 

i^M-t) - H{t)Mt) = (Hit) - i/o)V'oe-'^°*. (5) 
ot 



2 



Although Eqs. [3] and [5] are formally equivalent when 
the time propagation is exact, they will behave differently 
when it is approximated. In Eq. [31 the propagator acts 
on both the bound and the tunneling components, while 
in Eq. [5] it acts on the tunneling component alone. 

If the Hamiltonian H{t) is approximated by H{t), 
where H{t) = Hit) — C{t), the Hamiltonian error ({t) 
will cause the propagated wavefunction to deviate from 
the true wavefunction. Writing the true tunneling com- 
ponent iptix, t) = 4>t{x, t)-\-t{x, t) as the sum of the prop- 
agated component plus an error term, then propagating 
(j)t using the approximate Hamiltonian, so that 

(z^-i/)(V'h(a:,t) + (/.t(x,i)) = 0, (6) 

causes the error term to obey an inhomogcncous 
Schrodinger equation 

*^e(t) - H{t)e{t) = at){M^, t) + cj^m, (7) 

where the Hamiltonian error acts upon both the 
bound and the tunneling components of the propagated 
wavefunction to serve as a source term for the error. If 
the tunneling component is small relative to the bound 
state component, the dominant source of error is the 
Hamiltonian error acting upon the bound state. Thus, 
a small error requires that an eigenstate of the true 
Hamiltonian H{t), be accurately described by the ap- 
proximate Hamiltonian H{t) used in the propagation. 

In contrast, Eq. [5]allows for two approximations to the 
Hamiltonian - the approximation Hb — Hq — C^b used to 
find 0s, and the approximation Ht{t) — H{t) — C,t{t) 
used to propagate </>(. In terms of these Hamiltonians, 
the wavefunction error now obeys 



H{t)et{t) = CBi'B{x,t)+CT{t)iJt, (8) 



so that the error in the bound state Hamiltonian oper- 
ates on the bound state, and the error in the propaga- 
tion Hamiltonian acts on the tunneling component. The 
advantage of the inhomogeneous TDSE, then, is the free- 
dom to choose a different (presumably better) approxi- 
mation to the Hamiltonian when finding the initial bound 
state than will subsequently be used in the propagation. 
Doing so will reduce or eliminate the major source of 
error in Eq. [7] at the cost of finding a single, static eigen- 
function of Hb, while leaving the cost of propagation 
essentially unchanged. 



Gauge Transformation 

The use of an inhomogeneous Schrodinger equation has 
a long history in analytical approaches to strong field ion- 
ization, in particular strong field S-matrix theory 0, 



Here, the molecular potential is often ignored completely 
after the moment of ionization, so that the tunneling 
component can be expanded in a basis of Volkov states 
(j)t = J dpCp \p — eA{t)) , where p is the conjugate mo- 
mentum of the electron in the laser field, e = — 1 is the 
charge of the electron, and 

Cp^-i / dt'e-*^(P-^^(*"))'''*" {p ~ eA{t)\ VL{t') |0b) e-*^"* 

(9) 

or, in differential form, 

i^§i-lip-eAit)f)Cp = {p\ VLit') e^^^*'. 

(10) 

An apparent gauge dependence now arises from the 
question of what gauge VL{t) — H{t) — iJo is to be cal- 
culated in. As the Volkov states are stationary in the 
velocity gauge, it may seem natural to calculate Vl in 
this gauge as well. However, comparisons with the nu- 
merical TDSE have shown [ITf better agreement when Vl 
is calculated in the length gauge. Although perfect gauge 
invariance cannot be expected when ignoring the molec- 
ular potential, such discrepancies have been seen as well 
in the ionization of negatively charged ions, where the 
strong field approximation is better justified. As gauge 
invariance is an essential feature of a physical theory, a 
brief digression is warranted. 

As written, Eqs. [3] and [5] are manifestly gauge invari- 
ant. Rewriting the right hand side of Eq. [5] as 

(H - Ho)^oe-'''"* ^ (»^ - H{t))^oe-'^"\ (11) 
Eq. [5] becomes 

^g-^Mt) - H{t)Mt) = fix, t) = (*^ - H{t))4,^e-'^°\ 

(12) 

so that Eq. [3] and both sides of Eq. [T^] have the form 

if{x,t)^{i^^~H{t))^{x,t), (13) 

where 1^3(2;, t) can be thought of as an action residual 
which arises when ^j{x, t) does not satisfy the Schrodinger 
equation perfectly. If we now perform the gauge trans- 
formation -4}{x,t) ■0(a;,i)e'^^^'*\ A{x,t) -> A{x,t) ~ 
VA{x, t), 4>{x, t) — >■ 4>{x, t) + ^A(a;, t), then transforming 
the action residual according to (p{x,t) — >■ ip{x,t)e^^^^'*^ 
preserves equality. In Eq. ^ ip{x,t) — is unaffected by 
the change of gauge, while in Eq. [TH fix, t) is affected in 
the same way by transforming both sides of the equation, 
so that the theory as a whole is gauge invariant. 

The apparent gauge dependence which is seen in strong 
field S-matrix theory arises from calculating the source 
term in one gauge and the tunneling component in an- 
other. Thus, it is equivalent to transforming one side 
of Eq. [5] but not the other. It is therefore not a true 
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gauge transformation which causes the observed results, 
but rather a subtle modification of the effective source 
term. If we reverse the one sided transformation, so that 
both sides are calculated in the same gauge, the effective 
source term wiU not be the field free bound state ijjb{t), 
but rather ?/'f,(i)e~'^'^^'*^. For the specific case of strong 
field S-matrix theory, the action residual from Eq. 1101 
becomes 



iEt,t' 



(14) 

so that the effective source term when calculated in 
the same gauge as the tunneling component becomes 



i>bit)e~ 



eA(t)-. 



The improved accuracy of the 



"length gauge" calculation can now be explained by not- 
ing that the effective source wavefunction 



eA{t)-. 



obeys 



(Ho + E{t) ■ x)Mx,t), 



(15) 



(16) 



so that it is a quasistatic eigenfunction of the field-on 
Hamiltonian rather than a static eigenfunction of the 
field-off Hamiltonian. By incorporating the dressing of 
the bound state by the laser field into the source term 
rather than the tunneling component, a correspondingly 
smaller fraction of the dynamics must be treated using 
the inexact propagation Hamiltonian, giving a simple ex- 
planation for the improved accuracy. 



NUMERICAL PROPAGATION 
Least Action Propagator 

The inhomogeneous form of the Schrodinger Equation 
may be propagated using methods similar to those used 
for the homogeneous form. For this paper, we prom- 
gated both forms using the least action propagator [l8| . 
For each timestep, we expanded both the bound and the 
tunneling components of the wavefunction with respect 
to a spacial basis Xii^) and a temporal basis Tn{t), so 

that (l){x,t) = (l)o{x,t) + (f>t{x,t) = J2i,n^i,riXiix)Tnit). 

The expansion coefficients Din were found by minimiz- 
ing the action accumulated by the wavefunction over the 
chosen timestep. By choosing finite bases in space and 
time, this can be accomplished by solving a linear system 
of equations for Dm- In order to reduce computational 
cost, our spatial basis was chosen to yield banded ma- 
trices, while our temporal basis was constructed to allow 
for eigenvector decomposition. 



For a time dependent Hamiltonian H{t) and tempo- 
ral interval [to , to + At] , the action accumulated over the 
timestep is minimized when 



(zO, 



Hi 







(17) 



for alH, n such that Din is a free parameter (some values 
of Din will be fixed by the choice of initial or bound- 
ary conditions). Here Oij — (xilXj) is the spatial over- 
lap matrix, Qij = (T„| ^ \Tm) is the temporal derivative 
matrix, and Hijnm = {XiTn \ H{t) \xjTm) is the (time de- 
pendent) Hamiltonian matrix. For a static Hamiltonian, 
Hijnm = HijUnm, whcrc Unm = {Tn\Tm) is the temporal 
overlap matrix and Hij = {xi \ H\xj) is the Hamiltonian 
matrix. 

Separating the bound and tunneling components as 
in Eq [5] yields a linear system with a nonzero right 
hand side. In action language, the laser field causes the 
bound state to accumulate nonzero action, which must 
be canceled by the evolution of the tunneling compo- 
nent. Writing the tunneling component as (j)t{x,t) = 
X^i n CinXi{x)Tn{t) and the field free bound state compo- 
nent as 4>o = J2i n BinXi{x)Tn{t), the least action equa- 
tion becomes 

{iOijQnm ^ijnjn^Cjjn — AUijnm^jTn : (tS) 

where AHijnm = (Xi^nl AHsit) \XjTm)- 

In the velocity gauge, the kinetic energy term of the 
Hamiltonian is given by H{t) = i(p-eyl(t))2, where A(t) 
is the time dependent vector potential, with matrix ele- 
ments A„„. = (r„| A{t) \T„,) and Al„-^ = (r„| A^{t) \T„,) 
for its square. The squared term is often neglected in 
numerical treatments, but is required here for gauge in- 
variance and for use with exterior complex scaling. The 
least action linear system is now 



[iOijQnm ( AijU'i 



ij^nm ^'^ijAnni 



(19) 



where Rin is the inhomogeneous source term, V^j = 
(Xil^lXj) is the matrix element of the gradient and 
Aij — {xi\A\xj) is the matrix element of the Lapla- 
cian. As discussed earlier, the gauge in which Rin is 
found may be varied, with a corresponding change in the 
effective bound state. For this paper, we calculated Rin 
in both the velocity gauge 

and the length gauge 

O^JUnmRJm = {x^{x)Tnit)\e-'^'■'^'= xMx,t)) . (21) 

The linear system constructed in this way requires a 
large number of coefficients to be found at every timestep 
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- for a spatial basis set of size and a temporal basis 
set of size Nt, N^- Nt coefficients must be found at every 
timestep. We chose our spatial and temporal basis sets so 
as to decompose this problem into the iterated solution of 
[Nt — 1) sparse, banded linear systems for coefficients 
apiece, so that every timestep required computational 
effort on the order of NxNt. 

Spatial basis 

Our spatial basis was constructed following [l^ to al- 
low for infinite range exterior complex scaling. For a 
one dimensional problem, the region — cxd < x < oo is 
divided into a number of finite elements, plus two semi- 
infinite end caps. Within each element, the wavefunction 
is described using a primitive basis of low order Legen- 
dre polynomials. In each end cap, the primitive basis set 
is given by Laguerre polynomials times declining expo- 
nentials, so that the wavefunction is required to be at 
X = ±oo. In contrast to [3], here the same order poly- 
nomials were used in each element and both end caps. 

We dealt with the problem of reflections from the 
borders of the propagation region by employing infinite 
range exterior complex scaling. Beyond some scaling 
radius Xecs, we analytically rotated a; — Xecs + {x — 
a^ecs)e~'^, so that a wavefunction of the form e'*^^ will go 
to zero as x goes to infinity. This was accomplished by 
multiplying O^j O^C^ A^j -> A^e^^, Vije-''^ 
for all elements beyond Xecs- Because exterior complex 
scaling causes a derivative discontinuity at Xecs , it is nec- 
essary to enforce an element boundary at Xecs- In con- 
trast to we include a number of finite elements in 
the scaled region in addition to the semi-infinite end caps. 
The inner boundary Xf, of the end caps was chosen such 
that Q-kl^b-^o'iA ^ where k — \J2Emax, was less than some 
tolerance. 

Aside from enforced element boundaries at x = and 
X = zbxecs, the size of an individual element was de- 
termined using a WKB approach. Choosing an energy 
cutoff Emax and a maximum phase to be accumu- 
lated in an element, the size Ax of an element was given 
by A(/) = Ax^2{E„iax + \V{x)\ + \Fx\), where F is the 
maximum strength of the laser field. Constructing the 
spatial basis in this way ensures that a local Hamiltonian 
will yield sparse, banded spatial matrices, with a band- 
width determined by the order of the polynomials in the 
primitive basis set. Highly oscillatory functions can be 
treated by decreasing the size of the finite elements, or by 
increasing the order of the polynomials used inside each 
element. 

Wavefunction continuity across element boundaries 
was enforced by constructing linear combinations of the 
primitive basis sets. Within each element, we constructed 
left border functions BFL{y) = i(Po(y) — Pi{y)), right 
border functions BFbXu) = ^{Poiv) + Pi{y)), and inte- 



rior functions BFn{y) = (P„(y)-P„ mod 2(2/)) for n > 2, 
where y is a local coordinate equal to —I at the left 
boundary and 1 at the right boundary. Recalling that 
P/(-l) = (-1)' and Pi{l) = I, it can be seen that BFl 
is the only nonzero function at the left boundary and 
BFn is the only nonzero function on the right boundary. 
Wavefunction continuity was enforced by requiring the 
temporal coefficients Cm of the left border function in 
one element to equal the coefficients of the right border 
function in the element to the left, and vice versa. Border 
functions in the end caps were defined in an analogous 
way using Laguerre rather than Legendre polynomials. 

Temporal basis 

While the spatial basis was constructed to yield sparse, 
banded spatial matrices and to enforce wavefunction con- 
tinuity, the temporal basis set was constructed to reduce 
the size of the linear systems which must be solved at ev- 
ery timestep. Linear combinations of the primitive Leg- 
endre basis were constructed so that only one basis func- 
tion was nonzero at the beginning of the timestep, while 
the other basis functions had a tridiagonal overlap struc- 
ture. Using local coordinate r, equal to —I at the begin- 
ning of the timestep and -1-1 at the end, we set To{t) = 
i(Po(T)-Pi(T)) andT„(t) = i(P„(r)+P„_i(r)), where 
P„(t) is the Legendre polynomial of order n. Note 
that Toito) = 1, Mto + At) = 0, and r„(io) - 0, 
r„(to + At) = 1 for n > 1. Thus, the wavefunction 
at the beginning of the step is given by the coefficients 
of the zeroth order basis function, while the wavefunc- 
tion at the end of the pulse is given by the sum of the 
coefficients of all the nonzero basis functions. 

The size of the linear systems to be solved can be re- 
duced by performing an eigenvector decomposition of 
the residuals Pm, n > 1, which must be eliminated 
by the tunneling wavefunction. (Because Ciq are not 
free parameters, P^o will not in general be canceled by 
the minimum action solution.) Noting that for small 
timesteps A(t) and J^(t) are approximately constant, so 
that Anra ~ aUnm and A^j„j « bUnm, the residuals can 
be projected onto generalized eigenvectors of the Q and 
U matrices. Setting Qnm = Qnm and Unm = Unm for 
n,m> 1, U and Q have right and left eigenvectors 

Qw" = A"C/w" (22) 

and 

w^Q = X°'w°'U, (23) 

where A" is a generalized eigenvalue and w"^ Uv^ = 5ap ■ 
Ignoring (A„m - aUnm) and (^^„j - bUnm) for the mo- 
ment, if 5Cjyn = J2a^T'"mJ then 

{iO,jX''^i^Vl^iaV,j+iaVl^+bO,,+V,j))5Cf = ^<^P> 

n 

(24) 
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Because this method neglects the differences {Anm — 
aUnm) and {A\^ — bUnm), it wih yield an approxi- 
mate rather than exact solution for SC", so that ap- 
plying the correction Ci„ — Ci„ -I- SCin will not cancel 
the residuals completely. However, as a single timestep 
will typically be very short compared to the frequency 
of the driving laser, the variation of A{t) and ^^{t) 
over a single timestep will be very small, and repeated 
application of this procedure will yield rapid conver- 
gence to zero residual. We repeat this iteration until 
5ipt(x,t + At) = J2n^^'in has a norm less than some 
desired tolerance. After finding the minimum norm so- 
lution, we adjust the size of the next step so that the 
contribution of the last temporal basis function to the 
final wavefunction Sipn^a^ = J2i ^in^^^ Xi {x) has a norm 
less than a chosen accuracy goal. 

Having found the variationally optimum description of 
the wavefunction in the chosen basis, the timestep may 
now be adjusted to ensure that the wavefunction evolu- 
tion is well described by the temporal basis. We ensure 
this by choosing a timestep such that the last Legendre 
polynomial included in the primitive basis has a negligi- 
ble contribution to the total solution. If 

F^ = J2 ^™ (^"™- 1^") (25) 

n 

is the vector of spatial coefficients for the last Legendre 
polynomial and 

S^n^..=J2F*0,,F„ (26) 

i,j 

is its norm, then 

At' ^ accuracy goal ^„^^^_i , , 

gives the ratio of the new time step to the old. 

These design choices yield a variable stepsize, vari- 
ationally optimum propagator for a wavefunction de- 
scribed in a finite element basis. The order of the poly- 
nomials used to describe the wavefunction within each 
element, as well as the temporal order used to describe 
the wavefunction over a given timestep, can be chosen 
at will. The linear systems which must be solved are 
sparse, banded, and do not grow with increasing tem- 
poral order. For the calculations in this paper, we used 
polynomials of order 5 for both our spatial and temporal 
bases, and propagated with an accuracy goal of 10~^. As 
the spatial component (5'0n„ax of the last temporal basis 
function vanishes very rapidly with increasing rimax, the 
maximum step size was limited, not by the accuracy of 
the propagator, but rather by the accuracy of a polyno- 
mial expansion for the temporal evolution e*^"* of the 
source wavefunction. 



RESULTS 

We illustrate our treatment of strong field ionization by 
comparing the evolution of the tunneling component to 
that of the full wavefunction resulting from half a cycle of 
an intense laser pulse. For our model atom we use a one 
dimensional soft core potential V{x) = where 
Q — a — 1, whose ground state has energy Eg = —0.670. 
For laser electric field E{t) = Fq cos(wt), we set Fq = 0.05 
and oj = 0.0565 (equivalent to an intensity of 8.9 x 10^'^ 
W/cm^, consistent with many strong field experiments). 
The combination of ground state and laser pulse yielded 
a Keldysh parameter 7 — 1.31, placing it in the inter- 
mediate range between the tunneling and multiphoton 
limits. 

We compare the evolution of the full wavefunction to 
that of the tunneling component using three different po- 
tentials for the propagator. First, we verify that our ap- 
proach gives results identical to the homogeneous TDSE 
by comparing comparing their evolution using the exact 
propagator. We then demonstrate the effects of an imper- 
fect propagator by intentionally distorting the potential 
used in propagation. We perform one calculation using 
V(x) = (the well known strong field approximation) 
and one with V{x) = ~"y^^jf with a — > 2.0, so that 
the potential is asymptotically correct but differs from 
the correct potential at short range. Such a distortion is 
analogous to those which might arise in a single active 
electron mean field calculation, where the electron's in- 
teraction with the ion is correct at long range but differs 
from the true potential at short range. 

In Figures [T] and [21 we show that the inhomogeneous 
Schrodinger equation is identical to the homogeneous 
form for both the velocity and length gauge source terms. 
As discussed earlier, use of the "length gauge" source 
term is equivalent to using a quasistatic rather than a 
static wavefunction in the velocity gauge source term. In 
the present context, this means that the initial condition 
for the full wavefunction, when iptix, ti) = 0, must be the 
static bound state ^p{x,ti) = i/j^e"''^'*' for comparison 
with the velocity gauge calculation, and the quasistatic 
state 'ip{x,ti) = ^^e~*^9*'e~''^^^''^'^ for comparison with 
the length gauge. As seen in these figures, both source 
terms give results which are identical to the full wave- 
function in the region where ijjg{x) — 0. The apparent 
gauge dependence is thus reconciled: both gauges give 
results which are equivalent to the full TDSE, and differ 
from each other only due to the use of different source 
terms. 

The strong field approximation is shown in Figures [3] 
andlH Here setting V{x) — means that the propagator 
is incapable of treating a bound state correctly, so that 
the entire propagated wavefunction is projected into the 
continuum. As might be expected, propagating the full 
wavefunction gives results which bear no resemblance to 
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the true behavior. In contrast, the inhomogeneous ap- 
proach allows the propagator only to act on that por- 
tion of the wavefunction which differs from the initial 
bound state. Here the tunneling components have mag- 
nitudes similar to the true tunneling components, with 
peaks somewhat displaced due to the different potentials 
seen by the continuum wavefunctions. As the strong field 
approximation is commonly used in analytical treatments 
of strong field ionization, we show in Figure HJd the PPT 
tunneling component taken from 0. As in 0, only the 
exponential suppression due to tunneling is calculated, so 
that the ionization amplitude has only exponential accu- 
racy. The PPT wavepacket is much more localized than 
either the true tunneling component or the propagated 
strong field result, and has a considerably smaller ampli- 
tude. As interactions with the departing electron may 
cause dynamics within the parent ion, this may be an 
indication that more sophisticated treatments of strong 
field ionization will be necessary to describe processes 
such as multichannel ionization. 

For our third calculation, we observe the results of a 
distorted short range potential in Figures [5] and [5] In 
contrast to the strong field approximation, this propaga- 
tor treats the wavefunction correctly when it is far from 
the molecule, but incorrectly when the electron is nearby. 
Thus, propagating the full wavefunction will cause bound 
states to be projected into the continuum, but the con- 
tinuum wavepacket will evolve correctly for most of its 
journey. 

In the velocity gauge calculation, the full wavefunction 
shown in Figure [5^ has a very similar peak structure to 
the exact solution, but overestimates the amplitude of the 
wavefunction far from the origin. The tunneling compo- 
nent in Figure [SJj has amplitude comparable to the exact 
solution, but the maxima of the tunneling component are 
displaced somewhat from the exact solution. 

For the length gauge calculation, the incorrect poten- 
tial causes the initially quasistatic wavefunction to di- 
verge strongly from the exact solution in Figure EK, as 
a portion of the initial state is projected into the con- 
tinuum. In contrast. Figure [6]d shows a strong similarity 
between the exact and approximated tunneling compo- 
nents, with maxima which are once again displaced from 
the exact solution. 



CONCLUSIONS 

This paper has addressed the problem of strong field 
ionization from the perspective of the time dependent 
Schrodinger equation. Defining the tunneling component 
as that part of the wavefunction which differs from a 
precalculated bound state, we derived an inhomogeneous 
Schrodinger equation which governs its evolution. Our 
approach is manifestly gauge invariant, and we have ex- 
ploited this invariance to shed light on an apparent gauge 
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(b) 

FIG. 1. a) Tunneling component vs full wavefunction using 
velocity gauge source term b) Same calculation near the ori- 
gin. The tunneling component equals the full wavefunction in 
the region where the source wavefunction is zero, and differs 
where it is not. Lines show wavefunction magnitude, while 
color indicates phase. 



dependence which has long been puzzling in the context 
of the strong field approximation. 

While our approach is formally equivalent to the full 
time dependent Schrodinger equation in the limit that an 
exact propagator is used, it eliminates a major error term 
when the propagation is approximate. It thus addressed 
a common problem in the theory of molecular strong field 
dynamics, where the computational methods necessary 
to describe the bound state may be prohibitively expen- 
sive to be used in a time dependent problem. Use of 
the inhomogeneous Schrodinger equation may thus rep- 
resent a middle ground, allowing both the bound and 
the tunneling components to be treated using theoretical 
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FIG. 2. a) Tunneling component vs full wavefunction using 
length gauge source term b) Same calculation near the ori- 
gin. The tunneling component equals the full wavefunction 
in the region where the source wavefunction is zero, and differs 
where it is not. The initial condition for the full wavefunc- 
tion is the quasistatic eigenfunction rather than the field free 
eigenstate. Lines show wavefunction magnitude, while color 
indicates phase. 



FIG. 3. a) Full wavefunction propagated using exact propa- 
gator vs strong field propagator b) Velocity gauge tunneling 
component propagated using exact propagator vs strong field 
propagator. Whereas propagation of the full wavefunction 
fails completely for the strong field propagator, the tunneling 
component has comparable magnitude to that found using 
the exact propagator. Lines show wavefunction magnitude, 
while color indicates phase. 



approaches which are most appropriate to the task. 
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FIG. 4. a) Full wavefunction propagated using exact propa- 
gator vs strong field propagator for initially quasistatic state 
b) Length gauge tunneling component propagated using exact 
propagator vs strong field propagator and PPT wavepacket. 
Whereas propagation of the full wavefunction fails completely 
for the strong field propagator, the tunneling component has 
comparable magnitude to that found using the exact propaga- 
tor. The PPT wavepacket shows a markedly different magni- 
tude and structure from the full strong field calculation. Lines 
show wavefunction magnitude, while color indicates phase. 
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FIG. 5. a) Full wavefunction propagated using exact propa- 
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propagator vs propagator with modified short range poten- 
tial. Lines show wavefunction magnitude, while color indi- 
cates phase. 



(2010) . 

[15] L. Torlina, M. Ivanov, Z. Walters, and O. Smirnova, 

Physical Review A 86, 043409 (2012). 
[16] M. Spanner and S. Patchkovskii, Physical Review A 80, 

063411 (2009). 

[17] W. Becker and D. Milosevic, Laser physics 19, 1621 
(2009). 

[18] Z. Walters, Computer Physics Communications 182, 935 

(2011) . 

[19] A. Scrinzi, Physical Review A 81, 053845 (2010). 



9 




X 



(b) 



FIG. 6. a) Full wavefunction propagated using exact propa- 
gator vs propagator with modified short range potential for 
initially quasistatic state b) Length gauge tunneling compo- 
nent propagated using exact propagator vs propagator with 
modified short range potential. Lines show wavefunction mag- 
nitude, while color indicates phase. 



